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In this paper we compare different force fields that are widely used (Gromacs, Charmm-22/x- 
Plor, Charmm-27, Amber-f999, OPLS-AA) in biophysical simulations containing aqueous NaCl. 
We show that the uncertainties of the microscopic parameters of, in particular, sodium and, to a 
lesser extent, chloride translate into large differences in the computed radial-distribution functions. 
This uncertainty reflects the incomplete experimental knowledge of the structural properties of 
ionic aqueous solutions at finite molarity. We discuss possible implications on the computation of 
potential of mean force and effective potentials. 



I. INTRODUCTION 

The presence of water is characteristic to all biological 
systems. For any simulational study of biophysical phe- 
nomena, a proper description of water is a necessity. It is 
widely known that modelling water is a difficult task, and 
due to this many different models have been developed. 
In addition to water, the description of ions, in particu- 
lar Na + and CI - , is essential for biophysical systems. A 
simple example is the physiological liquid in the human 
body containing about 0.8 mol salt, i.e., of the order of 
one ion pair per 100 water molecules. 

Today, virtually all simulations use one of the water 
models either from the TIP or the SPC series, a recent 
review of water models is provided by Wallqvist and 
Mountain pj. For water, the influence of aspects such 
as density, treatment of long-range electrostatics and the 
choice of force- field have been studied intensively, see e.g. 
Ref. for a recent study. However, no systematic stud- 
ies, to the authors' knowledge, exist for ionic force fields. 
On one hand that is very surprising since various studies 
indicate that the force field can have a significant effect 
on the system properties even in the case of an implicit 
solvent pi]- In addition, it also known from experiments 
that the properties of ionic solutions may be significantly 
influenced even by small amounts of heavy isotopes of the 
ions 3. 

On the other hand, a practical obstacle for a system- 
atic comparison has been the evaluation of long-range 
electrostatics. In the past, there has been much work 
to include the effects of Coulombic interactions 0, IE S 
H, U E3, E3 without performing the computationally 
costly Ewald summations. Present day computational re- 
sources allow one to treat electrostatics properly |l3j and 
in, e.g., the case of lipid bilayers, this is a very important 
matter [l^. 

In this article we compare various commonly used 
force fields and a parametrisation for sodium chloride in 
combination with different commonly used water mod- 
els. For NaCl we used Gromacs, Charmm-22/x-plor, 
Charmm-27, Amber-1999 and OPLS-AA force fields, and 
the parametrisation by Smith and Dang |15| . The water 
models used for the aqueous solution were SPC, SPC/E, 



TIP3P and TIP4P. In previous studies some static prop- 
erties such as radial distribution functions have been 
studied, but in each case for one particular choice of force 
field only 0, H |1M GH 13 El In additiomthe effects 
of temperature [17| and salt concentration [18| have re- 
cently been studied. To our knowledge, there is no sys- 
tematic study regarding the effects of the force field and 
this study aims to fill that gap. 

When developing empirical force fields, one matches 
certain experimental quantities with their counterparts 
as determined by a numerical simulation. The choice is 
determined by the availability of high quality experimen- 
tal data and the physical significance of that quantity. 
The parameters of a force field thus depend critically on 
the choice of quantities that are being compared. Force 
fields are typically optimised for macroscopic parameters 
like the Gibbs free energy of hydration (i.e., placing an 
ion into a shell of water should give the same lowering of 
energy in the simulation as in an experiment) that are im- 
portant thermodynamic quantities and at the same time 
can be measured to a high accuracy. 

The thermodynamic properties are best complemented 
by structural properties. A natural way to quantify them 
is by using distribution functions of which the radial- 
distribution function gij(r) is the most commonly used, 
where r = \?i — fj\ stands for the separation between a 
particle of component i and of component j. It gives the 
probability of finding two particles at some distance r, 
taking account of density and geometric effects, and can 
be formally related to the potential of mean force be- 
tween two particles ^3- Aqueous NaCl is characterised 
by six different radial-distribution functions (since there 
are three components) but three of them cannot be de- 
termined experimentally at all, and of two only limited 
experimental information is available. We will return to 
this later in this paper. 

All force fields discussed in this paper have been pa- 
rameterised using the available experimental information 
on aqueous NaCl and therefore reproduce those experi- 
mental values reasonable well. In contrast, structural 
properties without known experimental values vary sig- 
nificantly between force fields. This uncertainty has be- 
come a significant problem recently since structural prop- 
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erties are essential in force field development for coarse- 
grained systems [2(1 0] . 



II. FORCE FIELDS 

In this article we compare different force fields. We 
restrict ourselves to force fields in the traditional sense, 
i.e., they have to be available in an electronic form and 
cover a wide range of systems. We specify the precise 
files used to obtain the parameters for our simulations. 
This information is relevant since often these parameters 
vary slightly between different sources even for the same 
force field. 

It should be noted that all of the force fields tested 
here were originally developed with the aim of describing 
proteins and nucleic acids. Description of ions is thus 
only a small part of their capabilities. For comparison, 
we have included one hand-optimised set of parameters 
for NaCl only HH- 

The different force fields and files are the following: 

Gromacs ("GROM"): Force field included in 
Version 3.1.3 of Gromacs j^- Available at 
http : / / www . gromacs . org/download/ index . php; 
file f f gmxnb . itp. The TIP4P water model for Gro- 
macs is available at http://www.gromacs.org/ 
pipermail/gmx-users/2001-November/ 
000152.html. For the systems discussed in 
this paper, i.e., water and NaCl, the Gromacs 
force field is identical to the Gromos-96 force 
field. Since Gromacs is the fastest MD program 
available, its default force field is used increasingly 
often. 

X-Plor/Charmm-22 ("XPLR"): Force field 
from x-plor distribution 3.851, available 
at http://atb.csb.yale.edu/xplor/; file 
parallh22x.pro. While this force-field is labelled 
as Charmm-22, the original Charmm-22 force 
field |2^| does not include ions, but they are in- 
cluded only in the x-plor distribution. X-plor [24| 
is one of the most versatile non-commercial pro- 
grams for protein simulations but is only able to 
use this force field. 

Charmm-27 ("CH27"): Available at 

https : / /rxsecure .umaryland . edu/research/ 
amackere/research . html; 

file par_all27_prot_lipid. inp. In comparison 
to Charmm-22, the more recent Charmm-27 [2^ 
contains parameters for ions in the files available 
at the website. Charmm-27 includes also other 
improvements for the description nucleic acids. For 
proteins, Charmm-22 is identical to Charmm-27. 
The Charmm-27 ion parameters are credited to 
Refs.ElandEll 

Amber-1999 ("AMBR"): The complete force field 
distribution for the Amber-1999 force field |28|| is 



available at http : //www . amber .ucsf . edu/ amber/ 
amber7 . f f parms . tar .gz. We used parameter 
file parm99 . dat and TIP4P water model from 
f rcmod.tip4p. New Amber-2002 force field in- 
cludes explicit polarisation terms, and thus falls 
outside the scope of this comparison. 

OPLS-AA ("OPLS"): The OPLS-AA force field H| 
is the only force-field in our list that is not part 
of a MD simulation package. Hence, there is no 
"official" file with the force field parameters. We 
chose the one included with Gromacs Version 
3.1.4, available at http://www.gromacs.org/ 
download/index.php in file f f oplsaanb . itp. 
One should note that other sources exist, e. g. 
http : //www . scripps . edu/brooks/ charmm_docs/ 
oplsaa-toppar . tar. 

Smith-1994 ("SMIT"): Hand-optimised set. Pub- 
lished in Ref. Il5l 

We performed simulations using the four standard 
water models, namely the rigid versions of SPC |3f)j |. 
SPC/E [a3, TIP3P [alEl and TIP4P 32]. For com- 
putational efficiency, we did not use the flexible versions 
of the water models. SPC/E and SPC differ only by 
the partial charges assigned to the atoms, so that the 
Lennard-Jones parameters are identical and thus need to 
be specified only once. 

The computer-readable files from the force field distri- 
butions typically contain one or more of the above water 
models. Whenever a water model was available in this 
way, we took the parameters from that file. Otherwise, 
standard parameters were used. This explains the (very) 
small differences than can be seen in Tab. |H between dif- 
ferent force field distributions for the same water model. 

The assignment of partial charges for the ions is triv- 
ial, and for the different water models it is well defined 
by the water model. The relevant parameters, since they 
are different for each force-field, thus are the ones describ- 
ing Lennard-Jones interactions. They can be specified in 
different ways, the two most common ones being 



V (r) = ^| - ^ 



4e 



(7)' 2 -(7 



(1) 



where the freedom of measuring energy in kcal or kJ re- 
mains. In addition, another common practise is not to 
specify all interaction parameters explicitly but to use 
the Lorcntz-Bcrthclot combination rules 



Cl + (72 



(2) 



where the indices 1 and 2 denote particles of type 1 and 
2, respectively. Table |I] lists the precise Lennard-Jones 
parameters used in our simulations. In addition, the ta- 
ble also indicates whether the parameter in question was 
specified directly by the force field or had to be computed 
via Eqs. {Q and/or @. 
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Gromacs 

Atom c 6 [kJnm 6 ] C12 [kjnm 12 ] e [^] a [A] 

CI 1.3804 ■ 10 2 1.0691 ■ 10 4 0.1064 4.4480 

Na 7.2059 ■ 10 5 2.1014 ■ 10 s 0.0148 2.5752 

O (S) 2.6171 ■ 10 3 2.6331 ■ 10 6 0.1553 3.1655 

O (3) 2.4889 ■ 10 3 2.4352 ■ 10~ 6 0.1519 3.1508 

O (4) 2.5543 ■ 10~ 3 2.5145 ■ 10~ 6 0.1549 3.1540 

CI— Na 9.9737 ■ 10 4 1.4989 ■ 10" 8 0.0396 3.3844 

CI— O (S) 6.0106 ■ 10 3 1.6778 ■ 10 5 0.1286 3.7524 

Na— O (S) 4.3426 ■ 10 4 2.3523 • 10 7 0.0479 2.8551 

CI— O (3) 5.8616 ■ 10 3 1.6135 ■ 10" 5 0.1272 3.7436 

Na— O (3) 4.2350 ■ 10 4 2.2622 • 10" 7 0.0473 2.8485 

CI— O (4) 6.4856 • 10" 3 1.9559 ■ 10~ 5 0.1284 3.8010 

Na— O (4) 4.4243 ■ 10~ 4 2.4446 ■ 10~ 7 0.0478 2.8646 



Charmm-27 



Atom 


c§ [kj nm 6 ] 


C12 [kJ 


nm 12 ] 


r kcal ] 


a [A] 


CI 


1.0999 ■ 


10~ 2 


4.8155 


10~ 5 


0.1500 


4.0447 


Na 


1.6169 ■ 


KT 4 


3.3284 


10" 8 


0.0469 


2.4299 


O (S) 


2.6171 


■ 10~ s 


2.6331 


• 10~ 6 


0.1553 


3.1655 


O (3) 


2.4912 ■ 


10" 3 


2.4364 


10" 6 


0.1521 


3.1506 


0(4) 


2.5543 


■ io- s 


2.5145 


• 10~ 6 


0.1549 


3.1540 


CI— Na 


1.6169 ■ 


10" 3 


1.8611 


IQ- 6 


0.0839 


3.2373 



CI— O (S) 5.6117 -lO' 3 1.2319 -lO' 5 0.1526 3.6051 
Na— O (S) 6.8542 ■ 10' 4 3.2868 ■ 10~ 7 0.0853 2.7977 
CI— O (3) 5.4847 ■ 10" 3 1.1892 • 10~ 5 0.1510 3.5976 
Na— O (3) 6.6750 ■ 10" 4 3.1500 ■ 10~ 7 0.0845 2.7903 
CI— O (4) 5.5514 ■ 10~ 3 1.2071 ■ 10~ 5 0.1524 3.5993 
Na— O (4) 6.7619 -IQ- 4 3.2028 ■ IP' 7 0.0852 2.7920 



OPLS-AA 



Atom 


C6 [kj nm 6 ] 




C12 [kJ 


nm 12 ] 


r kcal ] 


a [A] 


CI 
Na 
O (S) 
O (3) 
0(4) 


1.4654 


10" 


-2 


1.0886 


• 10" 


4 


0.1178 


4.4172 


6.3351 


10" 


-5 


8.6451 


■ 10" 


X 


0.0028 


3.3304 


2.6188 


10" 


3 


2.6352 


■ 10" 


(i 


0.1554 


3.1656 


2.4914 


10" 


3 


2.4367 


• 10" 


6 


0.1521 


3.1506 


2.5536 


10" 


3 


2.5121 


■ 10" 


(i 


0.1550 


3.1536 


CI— Na 


1.0227 


10" 


3 


3.4562 


• 10" 


(i 


0.0181 


3.8738 


CI— O (S) 


6.7301 


10" 


-3 


1.9991 


■ 10" 


5 


0.1353 


3.7914 


Na— O (S) 4.0810 


10" 


-4 


4.7915 


• 10" 


7 


0.0208 


3.2480 


CI— O (3) 


6.5798 


10" 


-3 


1.9314 


■ 10" 


5 


0.1339 


3.7839 


Na— O (3) 3.9820 


10" 


-4 


4.6110 


■ 10" 


7 


0.0205 


3.2405 


CI— O (4) 


6.6583 


10" 


-3 


1.9591 


• 10" 


5 


0.1351 


3.7854 


Na— O (4) 4.0311 


10" 


-4 


4.6810 


• 10" 


7 


0.0207 


3.2420 



X-Plor / Charmm-22 



Atom 


C6 [kj nm 6 ] 


C12 [kJ 


nm 12 ] 


e Si 


a [A] 


CI 


1.5362 • 


10~ 2 


9.3940 


10~ 5 


0.1500 


4.2763 


Na 


6.9284 • 


10~ 4 


2.8663 


lO" 7 


0.1000 


2.7297 


O(S) 


2.6171 


• 10~ 3 


2.6331 


■ w- e 


0.1553 


3.1655 


0(3) 


2.4913 • 


10" 3 


2.4366 


10" 6 


0.1521 


3.1506 


0(4) 


2.5543 


• 10~ s 


2.5145 


■ 10~ 6 


0.1549 


3.1540 


CI— Na 


3.7899 • 


10" 3 


7.0028 


10" 6 


0.1225 


3.5030 


CI— O (S) 


6.7840 


• 10- 9 


1.8004 


■ 10- 5 


0.1526 


3.7209 


Na— O (S) 1.3689 


■ 10~ :i 


8.9779 


■ 10- 7 


0.1246 


2.9476 


CI— O (3) 


6.6331 • 


10 -3 


1.7393 


10" 5 


0.1510 


3.7134 


Na— O (3) 1.3342 • 


10 -3 


8.6187 


10- 7 


0.1233 


2.9402 



CI— O (4) 6.7132 -lO' 3 1.7652 -lO' 5 0.1524 3.7151 
Na— O (4) 1.3513 -IP' 3 8.7593 -IP' 7 0.1245 2.9419 



Amber-1999 



Atom 


eg [kj nm 6 ] 




C12 [kJ 


nm 12 ] 


r kcal 1 
6 LmnlJ 


a [A] 


CI 


1.2170 ■ 


10" 


12 


8.8431 


10~ 5 


0.1000 


4.4010 


Na 


6.3072 ■ 


10" 


5 


8.5752 


10" 8 


0.0028 


3.3284 


O (S) 


2.6171 


■ 10 


-3 


2.6331 


• 10~ 6 


0.1553 


3.1655 


O (3) 


2.4904 ■ 


10" 


3 


2.4364 


10" 6 


0.1520 


3.1508 


0(4) 


2.5534 ■ 


10" 


3 


2.5116 


10" 6 


0.1550 


3.1536 


CI— Na 


9.2873 ■ 


10" 


4 


3.0945 


10~ 6 


0.0166 


3.8647 


CI— O (S) 


6.1201 


■ 10 


-3 


1.7946 


• 10~ s 


0.1246 3.7833 


Na— O (S) 4.0705 


■ 10 


-4 


4.7698 


■ lO' 7 


0.0207 3.2469 


CI— O (3) 


5.9839 ■ 


10" 


3 


1.7342 


10~ 5 


0.1233 


3.7759 


Na— O (3) 3.9722 ■ 


10" 


4 


4.5916 


10" 7 


0.0205 


3.2396 


CI— O (4) 


6.0564 • 


10" 


3 


1.7592 


10 -5 


0.1245 


3.7773 


Na— O (4) 4.0218 • 


10" 


4 


4.6612 


10- 7 


0.0207 


3.2410 



Smith-1994 



Atom 


C6 [kj nm 6 ] 


C12 [kJ 


nm 12 ] 


r kcal 1 
6 LmnlJ 


a [A] 


CI 


1.5790 ■ 


10" 2 


1.1458 


lO" 4 


0.1299 


4.4000 


Na 


2.8228 ■ 


KT 4 


4.7543 


10" 8 


0.1001 


2.3500 


(S) 


2.6172 ■ 


10" 3 


2.6338 


10" 6 


0.1553 


3.1656 


(3) 


2.4889 


■ lO' 3 


2.4352 


■ 10~ 6 


0.1519 


3.1508 


0(4) 


2.5543 


■ 10- 3 


2.5145 


■ 10~ 6 


0.1549 


3.1540 


CI— Na 


2.8223 ■ 


10" 3 


4.1711 


10" 6 


0.1140 


3.3750 


CI— (S) 


6.9705 • 


10" 3 


2.0424 


10" 5 


0.1420 


3.7828 


Na— O (S) 9.1847 ■ 


10" 4 


4.0406 


10- 7 


0.1247 


2.7578 



CI— O (3) 6.8132 -10 - 3 1.9730 ■ 10~ 5 0.1405 3.7754 

Na— O (3) 8.9383 -lO' 4 3.8693 -lO' 7 0.1233 2.7504 

CI— O (4) 6.8986 ■ 10' 3 2.0028 ■ 10' 5 0.1419 3.7770 

Na— O (4) 9.0590 -IP' 4 3.9352- IP' 7 0.1245 2.7520 



TABLE I: Parameters of the Lennard- Jones interactions for different force fields. The typeface of the numbers indicates where 
these numbers stem from. Boldface ("1.23") means that it is explicitly given by the force field in the specified notation. 
Underlined numbers (" 1.23 ") denote that the Lennard- Jones interaction parameters were given explicitly by the force field 
and a unit conversion (e.g. from kcal to kj) was necessary. Normal font ("1.23") means that the parameter in question was 
computed via the combination rule Eq. ||5J . Not all force fields specify all three water models. In case that one was missing we 
have taken the missing parameters (either directly or via the combination rule) from the Gromacs force field. This is indicated 
by italic font ("1.23"). Hydrogens do not participate in Lennard- Jones interaction, and the symbol after the O for oxygen 
stands for the water model ("S" for SPC and SPC/E, "3" for TIP3P, and "4" for TIP4P). 
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FIG. 1: Diffusion coefficients L>Na and Dei for sodium and chloride, respectively. Left: Labelling according to the ionic force 
field. Right: Labelling according to the water model used. 



In Sec. IIVI we present the results of our simulations, 
and show how the different force-fields differ in their de- 
scription of the NaCl properties. One important conclu- 
sion can, however, be drawn already from Tab. [I] The 
parameters for the different water models (SPC, TIP3P 
and TIP4P) differ only slightly, representing the current 
good knowledge of the properties of water. For (aque- 
ous) chloride, the differences are significantly larger, up 
to 10 % for the radius and up to 50 % for the depth of the 
attractive well of the Lennard-Jones interaction, reflect- 
ing the lack of high quality experimental input data. For 
(aqueous) sodium, there seems to be virtually no consen- 
sus on its properties. In the simulations, one can thus 
expect that the biggest differences will be in the Na-Na 
properties, followed by the Na-Cl interactions. 



slightly more than 10000 water molecules so that finite 
size effects are not expected. Lennard-Jones interaction 
was cut-off at 1 nm. The optimal choice for the cutoff 
length is not obvious and can vary between force fields 
(even between the ones for the ions and for water in the 
same simulation). For consistency, we decided to use the 
same cutoff in all simulations. For all of these systems, all 
relevant structures are on scales much smaller than 1 nm. 
Furthermore, all atoms are charged so that Lennard- 
Jones interactions quickly become negligible compared to 
electrostatic interactions, and the precise choice of cutoff 
does not matter as much as it does for other systems. 

The simulations described above presented a signifi- 
cant numerical task, and a total of approximately 25 000 
hours of cpu time was needed to complete them. 



III. SIMULATIONS 



IV. SIMULATION RESULTS 



For this study, we decided to include the three most 
commonly used thermostats, namely Berendsen |35j . 
Nose-Hoover |3(| |3?J and Langevin [33] . All of them are 
implemented into the Gromacs simulation software 
that was used for all of the computations presented in 
this paper. The target temperature was set to 298 K 
and particle-mesh Ewald (PME) was used for long-range 
electrostatics. The pressure was held constant at 1 bar 
using the Berendsen algorithm 35]. 

For each combination of ionic force-field, water model 
and thermostat a MD simulation was run. In addition, 
for each combination of water model and thermostat a 
reference simulation without ions was done. The total 
number of simulations added up to 84. A pre-production 
analysis showed that the systems needed slightly less 
than 0.5 ns to equilibrate. For each simulation run, we 
computed a 2 ns trajectory and only the second half of 
that was included into the analysis. 

The simulations were run at the physiological salt con- 
centration of 0.87 mol. The simulation box contained 



A. Dynamic properties 

The most common quantity to describe the dynamical 
behaviour of a system is its diffusion coefficient D. We 
have plotted the results for different forcefields and water 
models in Fig.^ The format for this plot, as well as the 
following ones, is the following: All results are plotted 
twice, using two graphs next to each other. The data 
points in both graphs are identical but in the left figure 
we have labelled them (i. e., picked symbols for the data 
points) according to the ionic force field whereas in the 
right figure we have labelled them according to the water 
model used. This way it is easy to see whether there is 
any systematic dependence on the ionic force field and/or 
the water model. 

The results for diffusion coefficients using Berendsen 
and Nose-Hoover were identical within statistical error, 
but using Langevin thermostat the diffusion coefficients 
were much smaller. Unlike the Berendsen and Nose- 
Hoover thermostats, the Langevin thermostat is not mo- 
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FIG. 2: Gibbs free energy of hydration for the different force fields per ion pair. 



mentum conserving, and thus we omit the results for the 
Langevin thermostat when computing the diffusion coef- 
ficients. 

The results of our simulations are depicted in Fig. ^ 
The experimental values at infinite dilution are given 
by D Na = 1.334 • 10" 5 cm 2 /s and D C i = 2.032 • 
10 -5 cm 2 /s |2^. From the figure it is seen immedi- 
ately that the dynamics is determined by the water model 
while the contribution of the ionic force field is negligible. 
The diffusion constant thus cannot be used to judge the 
quality of the different force fields. This is not surpris- 
ing since in aqueous systems the behaviour is dominated 
by the water molecules as they outnumber the ions by 
a factor of order 100. This is likely to be the case for 
other dynamic properties as well. A study of the dynamic 
properties of water models is beyond the scope and aim 
of this paper, and we refer to previous studies on this 
subject yi| and to the very informative webpage 0- To 
see the effect of the ionic force fields, we concentrate on 
static properties in the following. 

B. Gibbs free energy of hydration 

The most frequently used way to compute free energies 
is by the "slow-growth" method, also known as Kirk- 
wood's coupling parameter method |40l | or as thermody- 
namic integration. The presence of solute in a solvent is 
determined by the parameter A that modifies the Hamil- 
tonian H of the system. A = means that there is no 
interaction between solute and solvent whereas A = 1 
means that the normal Hamiltonian for the combined 
system is used. The change of Gibbs free energy upon 
introducing the solute into the solvent is then given by 



AG = 




where the averages are taken at constant pressure. (If 
the averages are taken at constant volume, this equation 
would yield the change of free energy.) The advantage 



of this approach is that the change of (Gibbs) free en- 
ergy is measured directly without the need to compute 
the (Gibbs) free energy of the solvent, thereby reducing 
the statistical error of the result. However, this method 
is computationally expensive since the integral has to be 
evaluated numerically by running simulations with dif- 
ferent values of A. While this is no significant problem 
when studying a single system, this is not feasible for a 
systematic study as presented in this paper. 

A different way to compute the free energy of a system 
is given by 




with U the potential energy of the system. This formula 
is easily understood by noting that it is conceptually 
identical to the computation of the chemical potential 
using insertion of a test particle, with the test particle 
being the entire system and the reference state being the 
vacuum. It offers the advantage that only a single simula- 
tion of the system with and without the solute is needed, 
not a large number of different simulations with gradu- 
ally increasing A. This formula is not as much used in 
the literature since for computing free energy differences 
it has to be applied twice to compute the free energies of 
the two comparison systems. The computed free energy 
difference then has a significantly larger statistical error, 
especially when only a small part of, for example, a large 
protein is mutated. 

We do not suffer from this problem as the number of 
ions is large enough to arrive at a statistically significant 
result. We thus apply Eq. (0J first to a simulation of 
aqueous NaCl and then to a simulation with the same 
number of water molecules (using the same water model 
and the same thermostat as in the first simulation) but 
without any ions. Since the two simulations are done 
at identical pressure and not at identical volume, the 
difference actually gives not the plain free energy but 
rather the Gibbs free energy G of hydration. 

The experimental values for the Gibbs free energy of 
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FIG. 3: Typical examples from our set of radial-distribution functions for Na-Na, Cl-Cl and Na-Cl (from left to right). All 
curves were computed using the Berendsen thermostat. 



hydration are —347 kJ/mol for chloride and —375 kJ/mol 
for sodium [3^, for a total of —722 kJ/mol. Older 
reported values are —363 kJ/mol for chloride and 
—406 kJ/mol for sodium It is assumed that iso- 

lated ions are to be hydrated, i.e., the crystal lattice of 
solid NaCl does no longer need to be broken up. (These 
two different concepts are frequently referred to as "hy- 
dration" versus "solvation".) It should be stressed that 
all those experimental values are at infinite dilution, i. e., 
all mutual interactions among the ions in the aqueous 
phase are ignored. 

The numerically computed Gibbs free energies G are 
depicted in Fig. |2 Since we used three different ther- 
mostats, this resulted in three different values for G for 
each combination of ionic forcefield and water model. We 
use those three different values to compute an error esti- 
mate that is marked by the error bar in the figure. Most 
of the error bars are that small that they hardly can be 
seen, confirming that application of Eq. is sufficient 
for this system. 

The experimental data given above for the Gibbs free 
energy of hydration is at infinite dilution and thus needs 
to be corrected to account for the finite molarity of phys- 
iological systems as treated in the simulations. We were 
unable to find direct experimental data for the enthalpy 
change of dilution but this quantity can be computed 
from the difference in the enthalpy of solution for the 
two different concentrations. For the parameters in ques- 
tion, the experimental values [42j lead to a correction of 
-2 kJ/mol. 

All numerically computed data in Fig. |3 are more neg- 
ative than the experimental data. This means that either 
the attractive mutual interaction of the ions is overesti- 
mated in the simulations, or that the shielding of that 
interaction due to the water is underestimated. It should 
be stressed that the parametrisation of all ionic force 
fields was done at infinite dilution, meaning that only 
a single ion was considered. The effect of finite concen- 
tration was thus not considered in the parametrisation 
process. A systematic numerical study of the concen- 
tration dependence of the Gibbs free energy of hydration 
could help to provide an understanding but no such study 
seems to have been published so far. Comparing experi- 



ment and the simulation results in Fig. fallows to judge 
the agreement between a force field and experiment - 
for precisely the molarity studied in this paper — but the 
picture might be completely different at other molarities. 



C. Radial-distribution functions 

Next, we compare the radial-distribution functions 
(rdf). There are three kinds of particles in the simula- 
tion, namely Na, CI and water, resulting in six different 
pairs and thus six different radial-distribution functions. 
For the purpose of computing the rdf's, the position of 
the oxygen atom is taken to represent the entire water 
molecule, and we will thus use the label "O" for water. 
For space reasons we will not discuss the water-water 
rdf's in the following. This is no relevant limitation since 
for the physiological concentrations they depend hardly 
on the ionic force field. (The complete set of rdf's can be 
found as supplementary material at our website |49).1 

In Fig. we present the rdf's for four different force 
field combinations. (For the complete set, see the sup- 
plementary material.) It is immediately obvious that the 
rdf's differ from each other in many aspects, such as the 
number of peaks, the relative and absolute heights of the 
peaks, and that those differences are significant. 

It would be impossible to present all rdf's by directly 
plotting them. To give a more systematic overview in 
a condensed way, we have computed the position and 
height of the first peak for all rdf's. For this, a Gaussian 
is fitted to the rdf in the neighbourhood of the peak. 
The results are depicted in Fig. 0] We first will discuss 
the simulation results and will then put them into the 
context of experimental results. 

From Fig.^Jit is seen that the Na-Na and Cl-Cl peaks 
are scattered widely (Na-Na being scattered more than 
Cl-Cl), and that there is no well-defined systematic ten- 
dency. This is the case for both the position and the 
height of the peak. Also for the Na-Cl results, the peaks 
are scattered widely but we want to point out one cu- 
riosity: While for the Na-Na and Cl-Cl peaks there is a 
dependence both on the ionic force field and the water 
model, the position of the Na-Cl peaks is independent 
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FIG. 4: Position and height of the first peak of the radial-distribution function for (from top to bottom) Na-Na, Cl-Cl, Na-Cl, 
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FIG. 5: Coordination number, computed by integrating the radial-distribution function between the water oxygen and the ion 
from zero up to its first minimum. 



of the chosen water model. An easy explanation could 
be that sodium and chloride ions are frequently in di- 
rect contact with each other so that the properties of the 
water model are of minor importance. 

The positions of the Na-0 and Cl-0 peaks differ less 
between force fields. (Please observe the different scale 
of the r-axis in the different subfigures.) The height of 
the peaks, in contrast, still varies by about a factor 3 
between the different force fields. 

We will now discuss these findings in the context of 
experimental limitations. Radial-distribution functions 
can, in principle, be measured by x-ray diffraction since 
they are related to the Fourier transformation of the 
structure factor. To get a strong enough signal, the sub- 
stances have to be of sufficiently high concentration. This 
makes the experimental determination of Na-Na, Na-Cl 
or Cl-Cl rdf's impossible. 

For Na-0 and Cl-0 rdf's, the signal is only weak at 
physiological concentrations. It is possible to determine 
the position of the peak of the rdf quite well since it is 
directly related to the wave length of the peak of the 
structure factor. The height of the peak of the rdf, how- 
ever, can be measured only with great difficulty. (We will 
return to this in the next section, when discussing the co- 
ordination number.) We should add that for sodium any 
kind of measurement is more difficult than for chloride 
since there exists only a single isotope of sodium, and 
consequently the difference method of isotope substitu- 
tion cannot be used. 

Determining suitable experimental values for the peak 
positions is aided by the experimental observation that 
they depend only very weakly on molarity of the salt. 
Surprisingly, the results for crystals and aqueous solu- 
tions are also very similar. Marcus 39] quotes values 
between 0.233 nm and 0.240 nm for the Na-0 peak posi- 
tion but earlier reviews include values up to 0.25 nm |43j . 
For chloride, the quoted values are between 0.30 nm and 
0.32 nm [33. Other values quoted in the literature are 
within these limits, except that also a Cl-0 peak position 



of up to 0.335 nm is reported [Zij . 

The simulation results in Fig. 0] are now easily under- 
stood. There are only two quantities that are susceptible 
to (decent-quality) measurement, namely the peak po- 
sitions of the Na-0 and Cl-0 rdf's. Consequently the 
predicted values for these two quantities vary only little 
between force fields. Taking the entire range of experi- 
mental values (i.e., not deciding by some criterion that 
one experimental value is "better" than the others), all 
force fields basically reproduce the experimental knowl- 
edge. 

The huge spread observed in the other peak positions 
and in all the peak heights in Fig. 0] is simply an echo of 
the lack of experimental data for these important struc- 
tural quantities, and no empirical force field can be better 
than the available experimental data. Any method that 
relies on mutual peak position of the ions therefore faces 
a problem, and no solution to this problem is in sight. 



D. Coordination number 

A different description of the water structure around 
the ions is given by the coordination number. Technically 
this quantity is computed by first determining the loca- 
tion of the first minimum of the radial-distribution func- 
tion and then integration the radial-distribution function 
from zero up to this point, 

"coord = -r / drr 2 g{r) , (5) 

►'water 7 

where VWater is the volume of one water molecule. The 
geometric meaning of n coor( j is that it gives the number of 
solvent molecules in the first hydration shell around the 
ion. It should be noticed that the coordination number 
is a purely geometrical quantity, and it does not imply 
that this number of water molecules is influenced in some 
way by the ion. 
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FIG. 6: Apparent molar volumes of the ions, and the change of the number of hydrogen bonds per ion pair. (The latter quantity 
is negative, meaning that hydrogen bonds are destroyed by the ions.) 



Experimental data on the coordination number is 
scarce. The reason for this is twofold: First, the am- 
plitude of the radial-distribution function is difficult to 
measure. The amplitude of g(r) cannot directly be in- 
ferred from the structure factor at the corresponding 
wave length but only follows indirectly from the normal- 
isation condition g(r) — 1 for r — > oo. Secondly, the 
position of the first minimum is less well-defined than 
the position of the first peak, and it is not always imme- 
diately obvious where to stop the integration. (This can 
already be seen from the examples of simulational rdf 's 
in Fig. El. 

The coordination numbers computed from our simula- 
tions are displayed in Fig. [5] First, the position of the 
first minimum was computed. For space reasons we de- 
cided against displaying these intermediate results and 
only show the resulting coordination numbers. (The po- 
sitions of the first minimum are between 0.305 nm and 
0.34 nm for Na-O, and between 0.375 nm and 0.44 nm 
for Cl-O.) 

The reported experimental coordination numbers for 
sodium are between 4 and 8 [43j. For chloride, the 
same source gives values for chloride of around 6 when 
measured for NaCl solutions but there are only few of 
those measurements. The coordination number of chlo- 
ride should change only slightly when the sodium is re- 
placed by some other cation. This should allow us to use 
also data for other salts, thereby getting better statistics. 
However, if this is done coordination numbers as high as 
II can be included [44), Given the large spread of the 
experimental values, the coordination number is not well 
suited to quantify the agreement between experiment and 
molecular dynamics simulations. 



E. Molar volumes 

Each ion occupies a certain volume determined by its 
radius. This quantity is of little practical relevance, how- 



ever, since the introduction of the ion changes the local 
structure of the surrounding water, and thus also the lo- 
cal density. These effects are included in the quantity 
known as molar volume. If the local compression of the 
water overcompensates the steric volume of the ion, the 
molar volume can become negative. 

The experimental molar volumes are —6.7 cm 3 /mol for 
sodium and 23.3 cm 3 /mol for chloride |39l l45l l4q . This 
gives a volume of 16.6 cm 3 /mol per ion pair. 

In a simulation using pressure coupling, the molar vol- 
ume per ion pair can directly be computed by determin- 
ing the average size of the simulation box, and from this 
subtracting the average size of the simulation box of a 
simulation of pure water. (It is essential to use the same 
water model and the same thermostat in both simula- 
tions.) This gives the molar volume per ion pair de- 
picted in Fig. Comparison with the experimental val- 
ues shows that Charmm-27, OPLS and Amber are per- 
forming a little bit less well in this test as do the other 
force fields. 



F. Hydrogen bonding 

Most ions are known to break the hydrogen bond net- 
work of the water around them. This can be explained 
by the water molecules aligning themselves for better in- 
teraction with the ion, at the expense of breaking the 
hydrogen bonds toward other water molecules. (On av- 
erage, a water molecule has about 1.55 hydrogen bonds 
in bulk water.) Experimental values for AHB, the change 
in the number of hydrogen bonds due to one ion, is —0.03 
for sodium and —0.61 for chloride [39ll47j. 

We have computed the change in the number of hy- 
drogen bonds by comparing two simulations, one with 
and one without the ions. Classical MD can only give an 
estimate of the number of hydrogen bonds, by counting 
the number of potential acceptors and donors that fulfil 
certain geometric conditions. The result is depicted in 
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FIG. 8: Life time of a contact between a cation and an anion, as well as the mean number of contacts that each ion has. 



Fig. It is immediately seen that all forcenelds yield 
values for AHB that are much more negative than the 
experimental values. 

An easy explanation can be offered for this observa- 
tion. All force fields discussed in this paper do not in- 
clude explicit energy terms for hydrogen bonds. Rather, 
hydrogen bond interaction is integrated out and is in- 
cluded in the partial charges inside the water molecule. 
This approach works very well if the water molecules are 
arranged in the same way as in bulk water. In any other 
arrangement, e. g. if ions are introduced, they are now 
easier to rotate as no hydrogen bond energy terms need 
to be broken. 



G. Cluster analysis 

Next we discuss the physical background behind the 
peaks in the ionic rdf 's in Fig. 0] For the Na-Cl inter- 
action a large peak is expected since the two ions are of 
opposite charge. For the Na-Na and Cl-Cl rdf's, this 
explanation does not hold. Especially if the peak is sig- 



nificantly higher than 1, this means that ions "like" to 
be at a certain distance from each other much more than 
to be away from each other as much as possible — even 
though they are strongly repelling each other through 
electrostatic forces. In earlier simulations of aqueous 
ionic systems pairing of chloride ions was found but it 
was later realised that this pairing is an artifact that 
disappears if long-range electrostatic forces are treated 
properly [110. 

Such problems can be ruled out here but there still 
is the question whether clusters of ions exist. In such a 
cluster, several cations (anions) can be close to each other 
because their mutual repulsion is compensated by the si- 
multaneous presence of several anions (cations) . For this 
reason we have performed a cluster analysis. We define 
a cluster as the set of all ions that are connected by dis- 
tances of 0.35 nm or less. From the radial-distribution 
functions in Fig. 0] it can be seen that 0.35 nm is a rea- 
sonable value, and that the precise choice does not effect 
the results. 

After doing this, we collected statistics on the number 
of ions in each cluster, computing the ratio p(N) of ions 
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that belong to a cluster consisting of N ions. In Fig.[7|we 
have plotted p(pair) = p(2) and p(cluster) = X^3P( n )- 
It can be seen that the different forcefields result in a 
large span for the results, and the probability for clusters 
of three or more particles spans more than two orders of 
magnitude. 

The different tendency for cluster formation is reflected 
also in the mean number of contacts that each ion has. 
This quantity is depicted in Fig. [S] Also depicted in this 
figure is the mean life time of each contact. This shows 
that if the life time is small, the mean number of contacts 
is also small but a large life time does not necessarily 
imply a large number of contacts. This can be explained 
by the interplay between the time scale for breaking a 
contact (i. c., the life time) and the time scale needed by 
two partners to find each other. 



H. Implications for the effective potential 

The radial-distribution function can be used to define 
different potentials. The most common one is the poten- 
tial of mean force Vpmf, defined by 

g(r) =cxp[-/?y P MF(r)] . (6) 

Although not immediately visible in the formula, the 
potential of mean force includes the direct interaction 
between two particles at fixed positions, and addition- 
ally the contribution from having a third particle at a 
fixed position provided particles 1 and 2 are already 
fixed [19j. In other words, the potential of mean force 
includes first order corrections to the pure pairwise po- 
tential. Differences in g{r) directly translate into differ- 
ences in VpmfM- 

If higher order corrections are included, a different 
kind of potential is found, termed effective potential 
V c s{r) 12] . The effective potential is defined by the con- 
dition that in the canonical ensemble it yields the de- 
sired radial-distribution function. Under the restriction 
that the Hamiltonian can be written as a sum of two- 
particle terms, the effective potential always exists and 
is uniquely defined (up to a physically irrelevant offset). 



If the full set of rdf 's is used for this process, the original 
microscopic potential is recovered. 

The value of the effective potential lies in that it allows 
one to integrate out degrees of freedoms. This process is 
referred to as coarse-graining. Of particular interest is 
to integrate out the solvent since this decreases the num- 
ber of particles and thus the computational burden by 
approximately a factor 100. Aim is thus to find three 
effective potentials, namely for Na-Na, Na-Cl and Cl-Cl 
interaction, that, when used in a MC or MD simulation 
without a solvent, reproduce the microscopic rdf's for 
Na-Na, Na-Cl and Cl-Cl computed including the sol- 
vent. 

As initial guess for V e g(r), the potential of mean force 
from Eq. (JSJ) is used. From the guessed effective poten- 
tials, the rdf's are computed in an MC simulation. If the 
rdf for particles of type i and j at distance r is larger 
(smaller) than the microscopic target rdf, the guess for 
V^\r) is decreased (increased). This process is re- 
peated until convergence is achieved. More sophisticated 
methods that use the four-particle correlation function 
to compute an improved guess for the effective poten- 
tial have been developed under the name inverse Monte- 
Carlo jl2j but they are not numerically stable enough for 
this problem. 

The effective potential for the rdf's from Fig. [3] are 
shown in Fig.[5| It is instructive to compare these two fig- 
ures since it demonstrates the difference between the po- 
tential of mean force and the effective potential. The Cl- 
Cl radial-distribution functions of both Gromacs+SPC 
and X-plor+TIP3P have a peak. According to Eq. © 
this translates into a dip in Vpmf- For V e ff, this dip is 
almost invisible for X-plor+TIP3P. This means that the 
peak in g(r) is (almost) solely due to interaction with 
particles that were not integrated out. In other words: 
Two chloride ions have a tendency to be close to each 
other due to attractive interaction with the same sodium 
ion. For Gromacs+SPC, Vpmf contains a big dip for the 
Cl-Cl interaction. This means that a big contribution to 
the peak in g(r) comes from interactions that have been 
integrated out, namely the water molecules. 

The important conclusion, however, is that not only 
the structure of the ionic solution, as given by g(r), de- 
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pends heavily on the force field used. Also the effective 
potentials show a large variation on the precise force field 
parameters. 

V. CONCLUSIONS 

In this paper, we have compared different force fields 
describing aqueous NaCl. To this end, we have computed 
several quantities, some of them of a thermodynamic ori- 
gin but being of structural nature. We can divide the 
computed quantities into two groups. In the first group 
are the quantities that have an experimental counterpart 
that is easily measurable. All force fields reproduce most 
of those values sufficiently well. In the second group are 
the quantities for which there exists no or very limited ex- 
perimental data. The different force fields predict vastly 
different values for those quantities. 

The Gibbs free energy of hydration is the most basic 
thermodynamic quantity for describing aqueous ionic so- 
lutions. It is a "purely thermodynamic" quantity if and 
only if the system is in the dilute limit. For finite con- 
centration of ions, the mutual interaction of ions becomes 
important, and thus the structural properties of the ions 
enter. Comparison of experimental and our simulational 
values for the Gibbs free energy of hydration shows that 
the difference between those values depends on the force 
field used, so in principle it can be used to judge the suit- 
ability of a force field for application at physiological salt 
concentrations. While the Gibbs free energy of hydration 



is an important fundamental thermodynamic property, 
its importance in most molecular dynamics simulations 
might be limited, however, since the ions are usually kept 
hydrated during the entire simulation. 

Structural properties are often more important, best 
described by the radial-distribution function. The posi- 
tion of the first peak of the radial-distribution functions 
for sodium-water and chloride-water is characteristic for 
the first hydration shell of those ions. All force fields re- 
produce the experimental values sufficiently well. This 
means that all force fields are able to describe well those 
physical effects that are governed by the hydration prop- 
erties in the immediate neighbourhood of the ions. 

For the ion-ion radial-distribution functions it was 
shown that the different force fields lead to significantly 
different results. Since there exists no experimental data 
for those properties, this does not imply that a few (or all) 
force fields were developed improperly but simply reflects 
the limited knowledge of such properties. The predictive 
power of any model for aqueous NaCl thus is very lim- 
ited if structural properties of mutual ion interaction are 
important. In particular this means that greatest care 
should be taken if the radial-distribution functions are 
used as input for some other application. 
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